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Abstract 

A new computational method is presented to resolve hydrodynamic interactions acting on solid 
particles immersed in incompressible host fluids. In this method, boundaries between solid par- 
ticles and host fluids are replaced with a continuous interface by assuming a smoothed profile. 
This enabled us to calculate hydrodynamic interactions both efficiently and accurately, without 
neglecting many-body interactions. The validity of the method was tested by calculating the drag 
force acting on a single cylindrical rod moving in an incompressible Newtonian fluid. This method 
was then applied in order to simulate sedimentation process of colloidal dispersions. 
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I. INTRODUCTION 



There are a number of useful systems consisting of small solid particles dispersed in host 
fluids. Among them, colloidal dispersions are most common to our daily life and are of great 
importance, particularly in the fields of engineering and biology . Colloidal dispersions 
have been reported to exhibit several unusual phenomena, such as long-range correlations 
in sedimenting particles 3], long-range anisotropic interactions in liquid crystal colloidal 
dispersions 0, transient gel formations during phase separations of colloidal su^ensions 



and electro-rheological effects in particle suspensions of nonconductive fluids 

Since the dynamics of colloidal dispersions are very complicated, it is extremely difficult to 
investigate their dynamic properties by means of analytical methods alone. Computational 
approaches are necessary in order to elucidate the true mechanisms of dynamic phenomena 
in a variety of situations. Colloidal dispersions, however, have a typical multi-scale problem. 
The molecules comprising host fluids are much smaller and move much faster than colloidal 
particles. From a computational point of view, performing fully microscopic molecular 
simulations for this kind of multi-scale system is extremely inefficient. An alternative, which 
is generally considered much better than microscopic simulations, is to treat host fluids as 
coarse-grained continuum media. 

Several numerical methods have been developed in an effort to simulate colloidal dis- 
persions. Two of the most well known methods are the Stokesian dynamics [7] and the 
Eulerian-Lagrangian method. The former is thought to be the most efficient Method, ca- 
pable of treating hydrodynamic interactions properly. Furthermore, it can be implemented 
as 0{Np) scheme for Np particles by utilizing the fast multipole method ^. However, it is 
extremely difficult to deal with dense dispersions and dispersions consisting of non-spherical 
particles by means of Stokesian dynamics due to the complicated mathematical structures 
used in Stokesian dynamics. On the other hand, the Eulerian-Lagrangian method is a very 
natural and sensible approach to simulate solid particles with arbitrary shapes. A number of 
kinds of tailor-made mesh, including unstructured mesh, overset mesh, and boundary-fitted 
coordinates, have been applied to specific problems, so that the shapes of the particles are 
properly expressed in the discrete mesh-space. Thus, in principle it is possible to apply this 
method to dispersions consisting of many particles with any shape. However, a numerical 
inefficiency arises from the following: i) re-constructions of the irregular mesh are necessary 
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at every simulation step according to the temporal particle position, and ii) the Navier- 
Stokes equation must be solved with boundary conditions imposed on the surfaces of all 
colloidal particles. The computational demands thus are enormous for systems involving 
many particles, even if the shapes are all spherical. 

Thus, our goal is to develop an efficient simulation method that can be applied to particle 
dispersions in complex fluids. Since host fluids are considered incompressible in such systems, 
an efficient simulation must address how to efficiently and accurately evaluate hydrodynamic 
interactions. As a first step towards this goal, we attempted to develop a method to simulate 
colloidal dispersions in simple Newtonian fluids. The reliability of this method was tested 
by calculating the drag force acting on a cylindrical object in a flow. Its performance was 
subsequently demonstrated by simulating the sedimentation processes of colloidal particles 
in a Newtonian fluid within a small Reynolds number regime. 

II. SIMULATION METHOD 

In order to overcome the problems arising at the solid-fluid interface in the Eulerian- 
Lagrangian method, rather than the original discontinuous rectangle proflle (interfacial 
thickness, ^ = 0) schematically depicted in Fig. ^ a smoothed proflle was introduced to 
the interface > 0). This simple modiflcation greatly beneflts the performance of numer- 
ical computations, compared to the original Eulerian-Lagrangian method for the following 
reasons. 

i) Regular Cartesian coordinates can be used for many particle systems with any particle 
shape, rather than boundary-fltted coordinates. The solid-fluid interface has a flnite volume 
(oc TTa'^^^^, with a and d as the particle radius and system dimension) supported by multiple 
grid points. Thus, the round particle shape can be treated in flxed Cartesian coordinates 
without difficulty. The simulation scheme is thus free from the mesh re-construction prob- 
lem that significantly suppresses the computational efficiency of the Eulerian-Lagrangian 
method. In addition, the simple Cartesian coordinate enables use of the periodic boundary 
conditions as well as the fast Fourier transformation (FFT). 

ii) At the interfaces, the velocity component in the direction normal to the interface of 
the host fluid must be equal to that of the particle. This kinetic condition is imposed in the 
Navier-Stokes equation, as the boundary value condition deflned for the solid-fluid interface 
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FIG. 1: An example of the smoothed profile (dashed line). The original rectangular profile is also 
shown for comparison (solid line). 



in typical methods. In our method, however, this condition is automatically satisfied by an 
incompressibility condition on the entire domain, which will be subsequently explained in 
detail. 

iii) The computational demands for this method include sensitivity to the number of grid 
points (volume of the total system), however it is insensitive to the number of particles. 
Thus, our method is thought to be suitable for simulating dense colloidal dispersions. 

The non-zero interfacial thickness ^ is the only approximation used in the present method. 
Thus, inter-particle hydrodynamic interactions can be fully resolved within the approxima- 
tion of the non-zero thickness in the present method. 

There have bee. .wo ..Ua. .ne.hods developed by d«e.e„t autho. SQ- The basic 
ideas of these methods is to use a fixed grid and represent the particles not as boundary 
conditions to the fluid, but by a body force or Lagrange multipliers in Navier-Stokes equa- 
tion. The essential difference in our approach to these two methods is to introduce a explicit 
diffuse interface in a smoothed profile. As long as a fixed grid is used, the moving boundary 
is inevitably represented as diffuse as grid spacing. The introduction of a explicit smoothed 
profile makes us to present a clear formulation of a numerical algorithm and is advantagous 
when it is applied to complex fluids [l], 12 1 

The lattice Boltzmann (LB) method 1^ has attracted much attention in recent years to 
simulate colloidal dispersions with hydrodynamic interactions ^|. The LB equation was 
proved to offer a faithful discretization of Navier-Stokes equation, and colloidal dispersions 
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are simulated in the Eulerian-Lagrangian manner. In practical viewpoints, the LB method is 
forumlated on a fixed Cartesian lattice and is well adapted to parallel computation. The LB 
approach, although the formulation is not intuitive and its treatment of moving solid-fiuid 
boundary is somewhat complicated, has several similar merits to the present method. 

The "fiuid particle dynamics" (FPD) method was proposed earlier and is similar to our 
method in spirit In this method, although a similar smoothed profile was adopted, 

there are several important differences between FPD and the present method. The most 
significant difference is that particles are modeled as a highly viscous fluid with viscosity rjc, 
much greater than the fluid viscosity r]s in FPD. This enables the rigidity of the particles to 
be sustained approximately by artiflcial diffusivity Ari(f){x, t) (Arj = rjc — ris ^ rjs) within the 
particle domain. While this model is physically correct, a practical problem remains in that 
a larger viscosity requires smaller time increments. In contrast, the present method treats 
colloidal particles as undeformable solids, {i.e., At] oo) thus no additional constraint 
arises in the numerical implementations. 

A. Basic working equations 

Colloidal dispersions are considered in a simple Newtonian liquid. The motion of the 
host fluid is governed by the Navier-Stokes equation with the incompressibility condition. 



where Uf is the fluid velocity, p is the fluid density. The stress tensor is represented by. 



where p is the pressure, rj is the fluid viscosity. 

The colloidal particles are assumed to be rigid and spherical, and their positions Ri are 
tracked in a Lagrangian reference frame. 



{dt + UfV)uf = -V-(Tf, 
V-Uf = 



(1) 
(2) 




(3) 



Ri 



(4) 



with the translational momentum equation 



(5) 
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and the angler momentum equation 



H 
i 1 



(6) 



where ii^, Vj, Oj, Mj, and Jj are the position, the translational velocity, the angular velocity, 
the mass, and the inertia tensor of the ith particle, respectively. The hydrodynamic force 
Ff and torque Nf acting on a particle can be obtained by integrating the stress tensor 
over the surface as 



where r is the relative position vector from the center of rotation to the colloid surface. 
Furthermore, F^^ is the force due to direct particle-particle interactions, and F^ — Mj(l — 
P*^)g is the buoyant force where p* = pc/ p is the mass density ratio of the particles to 
the host fluid, and g is the gravitational acceleration. Relevant dimensionless parameters 
in the above equations include the Reynolds number Re — UL/u, the Proude number 
Fr = U / \fgL-, the mass density ratio p*, and the volume fraction a. Here U and L represent 
typical velocity and length scales specific to the systems under consideration, respectively. 
The kinematic viscosity u = r]/ p and the mass density of colloidal particles pc are assumed 
to be constant. 

In typical methods, the above set of equations should be solved using proper boundary 
conditions defined at the solid-fiuid interface. In the present method, however, the solid-fiuid 
boundary condition is replaced with a body force and an incompressibility condition on a 
total velocity defined on the entire domain. 

B. Modified working equations 

Assuming a smoothed profile with a finite thickness ^ to the solid-fiuid interface, we here 
derive the body force which accurately takes the interactions between solids and fiuids due 
to the motions of colloids in an incompressible fiuid into consideration. The present study 
considers a mono-disperse system consisting of N spherical particles with radius a. The 




(7) 
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positions of the particles {-Ri, ■ ■ ■ , Rn} are first transformed to a continuous field 

N 

cf>{x,t) = Y^M^^t)^ (9) 

i=l 

using the ith particle's profile function (f)i{x) centered at R^. Several possible mathematical 
forms for (j)i{x) exist, however, some typical functions are listed in Appendix 1X1 

The continuum velocity field Up is defined for the solid particles using {Vi, ■ ■ ■ ,Vn}, 
{fii, ■ ■ ■ , SIat}, and as 

N 

(PUp{x, t) = {Vi{t) + X (33 - R,{t))} (Pi{x, t). (10) 

i=l 

The total (fiuid+particle) velocity field is then given by 

U{x,t) = {1 — (j))u f + (f)Up 

= Uf + (j){Up- Uf) . (11) 

Since the particle velocity field Up is constructed from the rigid motions of particles, V Up = 
is verified as 

V ■ Hp = V ■ ^ + X (a; - R,) ^1 

= Y,{y^ + ^^xi^- R^)} " V^, (12) 

i 

6i (V0i) - 0i (V0) 



02 

{V<Pi) (Pi - 0, (V0, 



0. (13) 



Assuming the incompressibility of the fiuid velocity Uj, the divergence of the total velocity 
is 

V ■ w = (V0) ■ {up - Uf). (14) 

The gradient of is proportional to the surface-normal vector and have a support on the 
interfacial domains. Therefore, the incompressibility condition on the total velocity V-n = 
means the solid-fiuid impermeability condition at the solid-fiuid interface. 

We are to derive the evolution of the total velocity u. To make the points clearer, we 
first consider the problem assuming that the motions of particles {Ri{t),Vi{t),fli{t)} are 



given. In Eq. only the fluid velocity Uf is to be solved. The evolution equation of the 
total velocity is splitted as 

{dt + u-V)u = -V ■ (15) 
dtu = 0/,, (16) 

where the stress tensor is 

a = -p*i + r]I^Vu+{VufY (17) 

By integrating Eqs (|15p with (jl7j) . the total velocity is predicted a.s u = u*. The pressure p* 
in Eq. (jl7p is determined to fulflU the incompressibility condition V ■ w = V ■ n* = 0. Then, 
the body force 0/^ is added to enforce Eq. (fTTj) and solid-fluid inpermeability condition. 
Therefore, the time-integrated body force 0/^ is determined as 

/ ds(j)fp = (j){up- u*) Vpph, (18) 

Jt P 

where the pressure Pp is determined to fulflll V ■ n = (V0) ■ {up — u*) = 0. By solving 

Eq. p6|) with the body force, Eq. ()18p. we flnally make Eq. (jlll) where the fluid velocity is 

(l-0)w/ = {1- ^)u* -^Vpph. (19) 

We note that the non-slip condition at the solid-fluid interface is fulflUed in this time evolu- 
tion of the total velocity. Since the viscous stress p7|) acts on the entire domain including 
the interfacial domain, the tangential velocity difference between Uf and Up is reduced. In 
other words, non-slip or slip condition can be imposed in the definition of the stress cr used 
in Eq. (jl5p . When compared to FPD their body force is (pfp = Ar^V-^ jVii + (Vw)^! 
with the artificial diffusivity Arj. In the limit of At] oo, the particle becomes rigid, how- 
ever this limit cannot be achieved by the numerical scheme used in FPD. In contrast to 
FPD, the body force 0/^ guarantees the rigidity of the solid particle without additional 
large artificial diffusivity. 

We complete the time evolution by deriving the hydrodynamic force acting on the 
particles. The hydrodynamic force is defined as the momentum fiux between the fiuid and 
solid. Thus, the hydrodynamic force is simply the counteraction from the fiuid. 



Jpdx. (20) 



In contrast to the force in Eq. (|7j) which is expressed as the surface integral, the above force 
is given as the volume integral. The volume integral is much advantageous in mesh-based 
discretization compared to the surface integral since no generation of body-fitted mesh is 
needed. 



C. Simulation procedure 

i) For a given particle configuration {R^}, velocity {V"}, and angular momentum {O"} 
{i = 1 . . .N), where the superscript n denotes the time step and h is the time increment, 
the fluid velocity at a time t = nh is predicted as 

dsV ■ y-cT - WW j , (21) 

under the incompressibility condition V ■ w* = which determines the intermediate pressure 

ii) At step n, the total velocity tt" should be equal to u.p within the particle domain 
and the surface-normal velocity components of particles and fluid should be match in the 
interfacial domain. Thus, it must be corrected by the body force defined by 

= {<^-{ul-u'))lh--^v,. (22) 

The correcting pressure is determined to make a resultant total velocity incompressible. 
This leads to the Poisson equation of Pp, 



Finally, we have 



VVp = pV-{0(w;-w*)}//i. (23) 

w" = w* + (24) 

= + Vv (25) 

iii) The hydrodynamic force and torque acting on each colloidal particle are now computed 
using the volume integrals, 

i^f = - I P<l>^fldx, (26) 

Nf = - J pix-R'Dx <jyj;dx. (27) 



and the velocity, angular velocity, and position of each colloidal particles at the current step 
n are updated to {n + 1) as 



Since the same 0/^ is used for both the host fluid and the colloidal particles fl^ . (fTTj) 
through the interface, no excess or shorts for solid- fluid interactions exist. The same type 
of solid-fluid interaction has been previously proposed in Ref. Q], though the treatment 
used in the present method is more general. It is important to note that the timing of the 
updates for fluid (j7T|) and particles (PHjl - lfHUI) has been shifted; the particles always go one 
step ahead of the fluid. This shift is primarily due to a technical issue related to the solid-fluid 
boundary condition. In general, the boundary value conditions, which is replaced with the 
force due to solid-fluid interactions in the present method, is necessary to update the fluid. 
Otherwise, in order to update both the fluid and the particles simultaneously, the implicit 
scheme for particles must be used. Consider the problem to integrate the (n — l)th step 
variables, {i?"^^, fi""^} with w""^, to nth step. The particle velocity of next step, 

u"^, in Eq. ^ is needed to update the total velocity u""'^ to before {J^^^ 
is updated. This situation requires the implicit treatment. The use of the implicit scheme 
complicates the algorithm and reduces the efficiency. The timing-shift is therefore necessary 
to realize the full explicit scheme described above. The initial condition, {-R", V", O"} with 
tt"~^, can be generically constructed. From the total velocity u^~^ satisfying the initial 
boundary condition given by {-R"~^, 0"~^}, the hydrodynamic force and torque is 

computed as surface integrals and the particle trajectory is integrated from (n — l)th to nth 
step. Together with the construction of the initial condition, the timing-shift algorithm is 
realized without loss of generality. 

In the above formulation, although the solid-fluid interaction force (|22j) was computed 
on the basis of the Euler scheme for simplicity of presentation, implementations of higher 
order schemes are straightforward for Eqs. (j28|) -(|3(J| ) . Furthermore, in order to update the 
host fluid in Eq. ()21|). no restrictions exist for the time discretization and any conventional 
scheme can be used for incompressible fluids. 




(28) 



(29) 



(30) 
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III. NUMERICAL RESULTS 



The present method has been apphed to two specific problems: The calculation of the 
drag force acting on an infinitely long cylindrical rod moving in a Newtonian fluid in order 
to check the validity. The method was also applied to simulations of many sedimenting 
particles in a two dimensional fluid in order to demonstrate the performance. In the present 
simulations, the Navier-Stokes equation was discretized with a de-aliased Fourier spectral 
scheme in space and a second order Runge-Kutta scheme (the Heun scheme) in time. For 
the colloidal particles, the velocity and the angular velocity were integrated with the Heun 
scheme, and the position was integrated with the Crank-Nicolson scheme. The external 
boundary condition on the edge of the systems was imposed in the same manner as the 
fluid-solid boundary condition on the particle surface. The simulation code is remarkably 
simple due to such unified treatment for all boundary conditions. 

A. Drag force on a cylindrical rod 

The drag force acting on an infinitely long cylindrical rod with radius a was computed 
by solving the Navier-Stokes equation around the rod in order to check the accuracy of the 
present method. Figure |2l shows a cross section of the geometry around the rod with finite 
thickness ^ at the interface. 

First, the effects of the finite thickness on the drag force are examined in the square 
domain of L^. An uniform stream U in x-direction was assigned to the edge of the domain 
as the boundary condition. Here the Reynolds number was defined by Re = 2aU / v. The 
drag coefficient was calculated as = F^/pU'^a, where the drag force was computed 
from Eq. (f^H|) for various values of A, a, L, and U. Figure El shows the relative error 
{CD{Re,^/A) - CniRe, ^/A = 0)} /CniRe, ^/A = 0) as a function of the interfacial thick- 
ness ^/A, where CoiRe,^/ A = 0) was estimated by extrapolating the measured curve of 
Co{Re,^/A) to ^/A — s> 0. The relative error in Co was observed to increase with increasing 
^/A, however, it tended to converge within 5% for several values of a/ A for ^/A = 1 and 
< Re < 20. Thus, ^/A = 1 was set for further simulations. 

Next, the drag coefficient Cd was calculated. The rod was fixed at the origin in the 
circular domain with radius L. The velocity at the external boundary r = L was set to 
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FIG. 2: Schematic representation of the cross section of geometry around the rod. A is the lattice 
spacing, a is the rod radius, and ^ is the interfacial thickness. The rod surface now has a finite 
volume ~ 27ra^ supported by several grid points on the fixed Cartesian coordinate. 

u{r = L,e) = U/{1 - 21og(a/L)) [{l - {a/Lf - 2\og{a/L)} - 2 {l - (a/L)^} cos ^e,] 
where and e,. are the unit vectors in x- and r- directions, respectively, and tan 6* = y/x. 
An analytical solution for the Stokes equation is known for this boundary condition, and 
the drag force is given by Fd = SnrjU /(I — 2 log(a/ Lj). The computed Cd using the present 
method as a function of Re is shown in Fig. |5 and is in good agreement with the theoretical 
Stokes law in Re < 1 within 5%. 

The accuracy of the present method using the finite interfacial thickness ^/A = 1 was 
determined to be acceptable for simulating colloidal dispersions for Re < 1 based on the 
numerical results above. 

B. Sedimentation 

The performance of the present method was examined by simulating sedimentation pro- 
cesses of mono-disperse particles in a two dimensional Newtonian fluid in a rectangular 
box surrounded by non-slip walls with p^, = 1.1 and a = 0.143 The dimensionless pa- 
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FIG. 3: The relative error in the drag coefficient Co as a function of the interfacial thickness ^/A. 

rameters were taken to be Re = 0.0916, Fr = 0.0512, where the settling velocity and 
the diameter of particle were taken as the characteristic velocity and length. Other com- 
putational parameters were chosen as A = 1, ,^/A = 1, a/A = 10, L^/A = 512, and 
Ly/A = 1024, where y-axis is in the direction of gravity. In order to prevent the particles 
from overlapping within the core radius ~ a, the force was added Ff ^ = -OEpp/dRi due 
to direct particle-particle interaction using the repulsive part of the Lennard- Jones potential 
Epp = 0AJ2f='i^ Ej=i+i [(2a/^i,)'' - i2a/R,jf] e{2"''a - Ri^), where ■ ■) is the step 
function and Rij = \Ri — Rj\. The direct interaction -Ff^ is not very important when 
the particles are moving around because the particles never overlap due to the lubrication 
effect, even without F[^. Fi gure El shows the lubrication force acting on two approaching 
rods computed using the present method. The lubrication force is always repulsive in this 
case, and thus prevents the rods from approaching each other. The strength of the repulsion 
increases with increasing velocity U. On the other hand, when the particles are stacked 
on the bottom wall during the later stage of sedimentation, F^^ is required to sustain the 
stacking against gravity. In fact, the repulsion vanishes for immobile pairs of rods. 

At the initial configuration, all the particles were placed near the upper wall and both 
fluid and particle velocities were set to zero, as depicted in Fig. IHfa). A typical snapshot 
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FIG. 4: Comparisons of the drag coefficient Cd (plus) from our method with the theoretical curve 
of the Stokes law (solid line). 

during sedimentation is shown in Fig. lEfb). Regions with swirled particles were observed, 
in which the particle velocities were highly correlated as a result of long-range inter-particle 
hydrodynamic interactions. A simulation with periodic boundary conditions in the horizon- 
tal (x-)direction was also conducted. In this simulation, swirls were still developed, however 
they were smaller than those observed with non-slip walls. The effect of confinement in 
the non-slip walls therefore enhances the velocity correlation. The computational demand 
required for the present simulation is less than one day of processing on a normal PC. 

IV. CONCLUDING REMARKS 

A new computational method has been developed to simulate particle dispersion in fluids. 
Utilizing a smoothed profile for solid-fluid boundaries, hydrodynamic interactions in many 
particle dispersions can be taken fully into account, both accurately and efficiently. In 
principle, the present method can be easily applied to systems consisting of many particles 
with any shape. The reliability of the method was examined by calculating the drag force 
acting on a cylindrical object in a flow. The performance of the method was demonstrated 
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FIG. 5: The normalized lubrication force acts on two approaching infinitely long cylindrical rods 
as a function of the nearest distance between the two surfaces. Different symbols denote different 



approaching velocities [/, ranging from 2.5 x 10 ^ to 9.6 x 10 ^ by C/ = 2.5 x 10 ^ x 2.61'' 



n 



... 11, which almost collapsed. The observed scaling behavior F = r]U f{h) with scaling function 
/(. . . ) is characteristic of Stokes flow, due to Reynolds numbers 2aU jv <\. 

to be satisfactory by simulating sedimentations of particles in a Newtonian fluid. 

Another primary benefit of using the smoothed profile arose when the method was ex- 
tended to colloidal dispersions in complex fluids with an internal degree of freedom, such as 
the molecular orientation or ionic density. In complex fluids, inter-particle interactions can 
be mediated by the internal degree of freedom of the fluid. In such cases, the fluid-particle 
interactions at the colloid surface could be more efficiently handled by utilizing a smoothed 
profile. Previous studies on particle dispersions in liquid crystal solvents demonstrate a strik- 
ing example of this efficiency 11, 3|- Although the hydrodynamic effects were neglected in 
these simulations, extensions to implement the hydrodynamic effects by incorporating the 
present method are currently underway. 
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FIG. 6: Snapshots of 240 colloidal disks sedimenting in a two-dimensional Newtonian fluid obtained 
using the present method. The magnitude of the host fluid velocity is indicated in color; change 
of color from blue to red corresponds to change of the fluid velocity from small to large. 

APPENDIX A: SELECTION OF SMOOTHED PROFILES 

The specific form of the smoothed profile should be selected according to the convenience 
of the physical modeling of systems under consideration. In the present study, an infinitely 
differentiable function with compact support was used. We adopted defined as 

Ux) = g{\x-R,\), (Al) 
aix) = h{{a + ^/2)-x) 

{expf-AVa;^) x > 0, 
X <0. 

where Ri, a, ^, and A were the position of the particle, the radius of the particle, the 
interfacial thickness, and lattice spacing, respectively. This choice is shown in Fig. |2l While 
this (p i^ciy appear somewhat complicated compared to other more simple choices, this 
(f) has the following benefits: i) three domains; solid, fluid, and interface are explicitly 
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separated, namely, = 1 is the solid domain {\x — Ri\ < a — ^/2), = is the fluid domain 
(a + ^/2 < \x — Ri\), and < < 1 is the interfacial domain {a — ^/2 < \x — Ri\ < a + ^/2), 
ii) high order derivatives of (pi with respect to x can be analytically calculated, and iii) due to 
its support-compactness, the integrals in Eqs. (j!26|) and (j27|) remain local, which contributes 
greatly to the efficiency of the computation. 
The second possible choice is 

^i{x) = - f tanh ^ ^ + 1 J . (A4) 

This choice was used in Refs jlll 12 , 1^ . This is also infinitely differentiable as well as 
analytically easy to handle. However, the support is not compact and the separation of the 
three domains is ambiguous. Furthermore, the ambiguity of domain separation tends to be 
more enhanced for higher order derivatives. For practical implementation, due to exponential 
decay of the hyperbolic function, a proper cutoff radius is adopted for the calculation of the 
integrals in Eqs. and (j2Z|)- 

The third possible choice is given by 

(j)i{x) = s{a - \x - Ri\), (A5) 
x< -e/2. 



Six] 



1 'sinf^ + 1) |a;| < ^/2, (A6) 



1 X > il2. 

which has the property of exact separation of the three domains, however, the second deriva- 
tive of is discontinuous at the fluid-interface boundary. Therefore, it is not recommended 
for computational models requiring derivatives higher than the second order of 0. 

The detailed choice of does not affect the results of the present simulations because 
only the first order derivative of is required in the present case. However, care must taken 
if higher order derivatives are required. 
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